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Abstract. - It is well-known that the original lattice Boltzmann (LB) equation deviates from 
the Navier-Stokes equations due to an unphysical velocity dependent viscosity. This unphysical 
dependency violates the Galilean invariance and limits the validation domain of the LB method 
to near incompressible flows. As previously shown, recovery of correct transport phenomena in 
kinetic equations depends on the higher hydrodynamic moments. In this Letter, we give specific 
criteria for recovery of various transport coefficients. The Galilean invariance of a general class of 
LB models is demonstrated via numerical experiments. 



Introduction. — The idea behind the discrete- 
velocity kinetic method is that the full details of the single- 
particle distribution function in the three-dimensional mi- 
croscopic velocity space is not entirely necessary in de- 
scribing the thermo-hydrodynamics of a fluid. Instead, 
much simpler kinetic systems in which the constituent par- 
ticles can only have velocities from a small discrete set can 
be devised to exhibit the same macroscopic behavior. This 
idea can be traced back to the works of Joule. It was later 
exploited by Broadwell [1,2] and more recently has lead 
to the development of the lattice gas cellular automaton 
(LGA) and lattice Boltzmann (LB) models [3-6]. The key 
question concerning the validity of this class of models is 
how closely the macroscopic thermo-hydrodynamics of the 
continuum kinetic system can be reproduced by the much 
simplified kinetic systems. 

It is well known that the original LGA [7, 8] has a non- 
Galilean-invariant hydrodynamics due to the artifacts in 
the advection term and the equation of state. The lattice 
BGK (LBGK) method [9,10] eliminated these artifacts by 
choosing a proper equilibrium distribution function in the 
single-relaxation-time (BGK) collision term. The Navier- 
Stokes hydrodynamics at the small-Mach-number limit is 
reproduced. Despite its phenomenal success in model- 
ing near-incompressible athermal flows [11-14], the LBGK 
method is still not fully Galilean invariant due to a "cu- 
bic" velocity dependency of the viscosity [15], which is one 
of the factors that limit the valid domain of LB method to 
near-incompressible flows. Several attempts [16-18] have 
been made to overcome this short-coming using more ve- 



locities and higher order terms in the Mach-number ex- 
pansion of the equilibrium distribution. The success is 
however limited for that not only the process of obtaining 
the high-order correction is tedious and model-dependent, 
but also the Galilean invariance is in some cases only par- 
tially restored. 

Recently, the LB method was re-formulated as a Her- 
mite series solution to the continuum BGK equation [19, 
20], in a way similar to the Grad 13-moment system. Un- 
der this formulation, the Chapman-Enskog [21] procedure 
can be applied to the LBGK system in the Hermite space, 
revealing that the convergence of LB to continuum BGK 
relies on the few leading moments of the equilibrium dis- 
tribution function. Provided that these moments agree 
with those of the Maxwell distribution, the macroscopic 
hydrodynamics of the LBGK equation is the same as that 
of the continuum BGK equation, which is known to be 
fully Galilean invariant. Once scrutinized in this frame- 
work, the commonly known LB models are insufficient in 
meeting the above conditions of retaining hydrodynamic 
moments, and therefore introduce errors such as the well- 
known cubic velocity dependency of the viscosity. This 
analysis is general enough so that it can be applied to 
the energy equation, the diffusion equation in a multi- 
component system, and higher hydrodynamic approxima- 
tions beyond the Navier-Stokes level to give convergence 
criteria of other transport coefficients at the Navier-Stokes 
level and beyond. 

In this Letter, we define systematically the necessary 
and sufficient conditions for the LBGK hydrodynamics 
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to converge to that of the continuum BGK. Essentially, 
these conditions are that sufficient hydrodynamic mo- 
ments must be retained in the equilibrium distribution 
of the LBGK system and accurately represented by the 
discrete velocities. We also present numerical simulation 
results that validate these conditions. In the next sec- 
tion, we first give a more elaborated derivation of the 
Chapman-Enskog calculation in Hermite space, and then 
point out the convergence condition as the direct conse- 
quence. These convergence conditions are then verified nu- 
merically by directly measuring the transport coefficients 
in LB simulations. Further discussions are offered in the 
last section. 

Convergence conditions. — In a previous publica- 
tion [20], the asymptotic convergence of LB to the con- 
tinuum BGK was briefly shown using Chapman-Enskog 
approximation procedure in Hermite space. Same as in 
continuum kinetic theory, the macroscopic equations of 
LB are essentially determined by the leading tensorial 
hydrodynamic moments. With discrete velocities, these 
moments are expressed as weighted sum of tensors con- 
structed from discrete velocities, which are not always 
isotropic as their continuum counterparts are. The insuf- 
ficient isotropy of certain moments has long been recog- 
nized as one of the reasons responsible for the deviations 
from the continuum hydrodynamics. By expanding the 
distribution function in Hermite polynomials, this tedious 
calculation can be tremendously simplified and extended 
to higher orders. 

We first note that the hydrodynamic equations are given 
by the conservations of mass, momentum and energy. 
Omitting external forces, they are: 

g+,V.«-0 (1) 

pf + V-P = (2) 

Pj f +P- Vw + iv-5 = 0, (3) 

where d/dt represents substantial derivative and the den- 
sity p, fluid velocity u, internal energy density per mass e, 
pressure tensor P and heat flux S are all velocity moments 
of/: 
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where / is the single-particle distribution function and £ 
the microscopic velocity. Clearly, two distribution func- 



tions with the same leading velocity moments will yield 
the same hydrodynamic equations. 

For closing Eqs. (l)-(3), the task of obtaining the ap- 
proximated distribution function using the few hydrody- 
namic moments was in the center of the development of the 
kinetic theory in the last century. In the Chapman-Enskog 
calculation, a successive sequence of approximated distri- 
bution functions are introduced in the order of Knudsen 
number [21]. For the following Boltzmann equation with 
the BGK collision model [22]: 
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where r is the relaxation time, the i-th approximation of 
the distribution is given by the recursive relation: 
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with /(°) is the Maxwell-Boltzmann distribution. On sub- 
stituting /W into Eqs. (l)-(3) for i = 0, 1 and 2 and using 
results of the previous approximation, the hydrodynamic 
equations at the Euler, the Navier-Stokes, and the Burnett 
levels are obtained [23]. This sequence of approximation 
quickly becomes formidably complex beyond the first two. 

Eq. (10) implies a simple recursive relation among the 
moments of the distribution functions at various approx- 
imation level. We note that the velocity moments of 
the distribution function are essentially its Hermite co- 
efficients [24]. Let be the n-th Hermite coefficient of 
the distribution /W. As previously shown, on substitut- 
ing the corresponding Hermite expansions into Eq. (10), 
the Hermite coefficients of the first approximation can be 
explicitly expressed using the Hermite coefficients of the 
zero-th approximation, namely: 
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^From the above equation, an important conclusion that 
can be immediately drawn about the BGK equation (9) 
is that the leading n moments in the first approximation, 
which gives the Navier-Stokes hydrodynamics, is given by 
the leading n+1 moments of J*- -*. The form of the hydro- 
dynamic equations is completely determined by the few 
leading moments of /(°). For instance, at the Navier- 

(2) 

Stokes level (i = 1), P is determined by a\ , which is 
in turn determined by the expansion coefficients of /(°) 

(3) 

up to Oq . Therefore, the momentum equation is correct 
as long as /(°) agrees with the Maxwellian in the first three 
terms of their Hermite expansions. 

For a finite sum of Hermite series, the leading moments 
are completely determined by the function values at a fi- 
nite set of abscissas. This property, together with the fact 
that the hydrodynamic equations are determined by the 
finite number of leading Hermite terms of f^°\ allows a 
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discrete-velocity kinetic system to have the same macro- 
scopic hydrodynamics as the continuum system when their 
leading moments are the same. Consider a Unite sum of 
Hermite series as /(°): 

/^^(flE/o'^^' ( 12 ) 

Let M be the highest order of the moments that deter- 
mines the hydrodynamic equations of interest, e.g., M = 3 
if the momentum equation at the Navier-Stokes level is of 
concern. The leading M (< N) moments are in the form 
of 

J w(€)p«)d£, (13) 

where is a polynomial of an order < M + N. This 
integral can be exactly evaluated using the values of the 
integrand on a finite set of velocities, {£, : i = 1, • • • , d}, 

as: 

^ d 

u(Z)p(Z)dt = J2 w iP(Zi)> ( 14 ) 
J »=i 

if and only if & are the abscissas of a Gauss-Hcrmite 
quadrature of a degree of precision Q > M + N, and Wi 
the corresponding wights. It can now be concluded that 
the sufficient conditions for the LBGK system to have the 
correct hydrodynamic equations are: 

1. The equilibrium distribution retains the necessary 
moments. 

2. The discrete velocity set allows the moments to be 
exactly evaluated using finite function values. 

Quantitatively, these two conditions are: 

N>M and M + N<Q. (15) 

Since Q is finite for any finite set of velocities and N < 
Q, N must be finite, i.e., the equilibrium distribution in 
LBGK is the sum of a finite Hermite series of an order 
which is smaller than the degree of the quadrature. The 
difference between Q and N is the maximum order of the 
moments whose dynamics will be the same as that in the 
continuum system. 

The requirements for the momentum equation to be 
fully Galilean invariant at Navier-Stokes level is imme- 
diately clear. Since M = 3 in this case, we have N > 3 
and Q > 6. Namely must retain all Hermite terms 
up to the third order and the discrete velocity set must 
form a 6th-order accurate Hermite quadrature. Neither of 
the original LBGK models [9,10] satisfies this condition as 
only the second order terms are retained and the velocity 
set is only 5-th order in both cases (Q = 5). As the veloc- 
ity set is sufficient to support the second-order terms, the 
error introduced in both models are due to the last term 
on the right-hand-side of Eq. (11), which manifests as the 
"cubic" velocity-dependence of the viscosity [15,20]. 



The analysis above can be applied to the case of energy 
equation where full recovery of the Navier-Stokes energy 
equation requires N > M > 4, and Q > 8. Using a 
third-order expansion (N = 3), or a velocity set that only 
supports third-order moments (6 < Q < 8) could yield 
otherwise correct thermal LBGK models with a velocity- 
dependent thermal diffusivity. 

It is worth pointing out that when conditions (15) are 
satisfied, the equilibrium function of Eq. (12) results in the 
exact hydrodynamics equations regardless of Mach num- 
ber. Hydrodynamic moments up to a given order can be 
shown to have the same dynamics of their continuum coun- 
terparts. These conditions are also necessary if a finite or- 
der polynomial is used as the the equilibrium distribution. 
Nevertheless, it was shown previously [25,26] that discrete 
velocity models can be constructed with exact conserva- 
tion laws and correct higher moments in the small Mach 
number limit. 

We note that for a given set of velocities and weights, 
Q = {(£i,Wi) : i = 1, ■ ■ ■ ,d}, the isotropy of the following 
tensors play a critical role in the analysis of LGA and 
LBGK [8]: 



d 




Here, we point out that Q corresponding to a Q-th order 
Gauss-Hermite quadrature is equivalent to the condition 
that: 

*•>-{*■> :™„ • ™ 

where 5^ is the rank-n isotropic tensor. Noticing that 
the above is always true in continuum, i.e., after defining 

Pn(£) = ' We haV6: 

n times 

the forward equivalence is trivial since Vn < Q, p n (£) is a 
polynomial of a degree < Q. Because of Eq. (14), 

= . (19) 

The backward equivalence is also straightforward. As 
■ n < Q} is a, linearly independent set and forms 
a complete basis of the Hilbert space of all polynomials 
of a degree < Q, any such polynomial can be written as 
a linear combination of p n (£)- By direct integration, it 
can be shown that Eq. (14) is true as long as Eq. (17) is 
true. Condition (17) is also proved by a slightly different 
procedure elsewhere [27]. 

Numerical simulations. — To numerically demon- 
strate the theory outlined in the previous section, here we 
present simulation results of some simple benchmarking 
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problems using LBGK models with equilibrium distribu- 
tions and velocity sets of different orders. All velocity 
sets are in three-dimensions, and they are the well-known 
D3Q19 model (Q = 5) [9], the seventh-order (Q = 7) 
39-velocity quadrature £ 3 19 7 of Ref. [20] (Q = 9) and 
a ninth-order 121-velocity quadrature -E39 1 . The two- 
dimensional projection of the last model was used in a 
previous study [28]. Here given in Table 1 are the full 
detail of the three-dimensional version. All simulations 
are performed on a 1 x 1 x 100 grid in Cartesian coordi- 
nates {x,y,z} with periodic boundary condition used in 
all directions. 

Table 1: The abscissas and weights of a ninth-order accurate 
Gauss-Hermite quadrature formula in three-dimensions. Here 
p is the number of abscissas in the symmetry class. The sub- 
script fs denotes a fully symmetric set of points. 





V 


W a 


(0,0,0) 


1 


0.03059162202948600642469 


(r, 0,0) FS 


G 


0.09851595103726339186467 


(±r, ±r, ±r) 


8 


0.02752500532563812386479 


(r, 2r,0) FS 


24 


0.00611102336683342432241 


(2r, 2r,0) FS 


12 


0.00042818359368108406618 


(3r,0,0) FS 


6 


0.00032474752708807381296 


(2r, 3r, 0) FS 


24 


0.00001431862411548029405 


(±2r, ±2r,±2r) 


8 


0.00018102175157637423759 


(r, 3r,0) FS 


24 


0.00010683400245939109491 


(±3r, ±3r, ±3r) 


8 


0.00000069287508963860285 
r = 1.19697977039307435897239 



To measure the velocity-dependence of the viscosity, we 
simulated the one-dimensional shear wave problem where 
both density and temperature To are constants initially. 
The initial velocity is given by: 

u = u + aoe x sm(2nz/L z ), (20) 

where Uo is a homogeneous constant translational velocity 
field, do the small initial amplitude of the shear, e x the 
unit vector in the x direction, and L z = 100 the periodic- 
ity in z direction. The initial distribution function / is as- 
signed with its equilibrium value based on the initial den- 
sity, temperature and velocity. When the Navier-Stokes 
equations are fully satisfied, the amplitude of the shear 
wave shall decay exponentially independent of Uq with a 
decay rate proportional to the viscosity. The viscosity is 
measured through the decay rate of the amplitude. 

Shown in Fig. 1 are the measured viscosity, normalized 
with respect to its theoretical value, as a function of the 
longitudinal and transverse components of Uq. Here the 
magnitude of Uq are expressed in terms of the correspond- 
ing Mach numbers. On the top are the results using the 
19 point (D3Q19) model with second and third order equi- 
librium distribution function (N = 2 and 3 respectively). 
The dependency on the translational velocity is significant 
in both cases except for small Mach numbers. For the case 
of N = 2, the error is caused by the missing third order 




Fig. 1: Velocity-dependence of viscosity of the fifth order 
D3Q19 model (left) and the seventh order 39-velocity model 
(right) using second- and third-order equilibrium distribution 
function. Here, vq is the theoretical value of viscosity at zero 
mean velocity and c = y/To is the isothermal sound speed. 



term in the equilibrium distribution function. For the case 
of N = 3, the 19- velocity quadrature is not sufficiently ac- 
curate for the third-order moments and results in an error 
in another form. On the bottom are results using the 39- 
velocity model. To be seen is that the velocity dependence 
of the viscosity is removed once the third order term is in- 
cluded in the equilibrium distribution function. 

The thermal diffusivity is measured in a similar fash- 
ion. Instead of velocity, a initial sinusoidal temperature 
perturbation is imposed while both the velocity and pres- 
sure are uniform. Fig. 2 shows the thermal diffusivity as 
functions of translational velocity uq and the background 
temperature To- Here uq is taken to have the equal com- 
ponents in three directions. As shown in Fig. 2a, with the 
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Fig. 2: Velocity-dependence of thermal diffusivity for the 121 
point model. The results in (a) and (b) are from the third-order 
and fourth-order equilibrium distribution function respectively. 
Here, ko is the theoretical value of diffusivity at zero mean 
velocity. 



Discussion. — In this Letter, we give the quantitative 
conditions for the LBGK system to have the same macro- 
scopic hydrodynamics as the continuum kinetic system de- 
scribed by the BGK equation. The degree of precision of 
the velocity set as a Gauss-Hermite quadrature together 
with the order of the Hermite terms retained in the distri- 
bution function determines the order of the hydrodynamic 
moments that will have the correct macroscopic behavior. 
The velocity dependency of the viscosity in the commonly 
used LBGK system is identified as an error due to both the 
insufficient truncation of the distribution function in mo- 
ment space and the insufficient isotropy of the 9-velocity 
model. Once sufficient moments are retained in the equi- 
librium distribution and accurate quadrature is used, the 
hydrodynamic behavior of the LBGK asymptotically ap- 
proaches to that of the BGK equation in continuum. Ki- 
netic models for thermal fluids, fluid mixtures and fluids 
beyond Navier-Stokes hydrodynamics can be constructed 
with guaranteed Galilean invariance. 

Centered in the effort of restoring full Galilean invari- 
ance in LBGK system is the search for velocity set that 
makes isotropic tensors. This effort is shown here to 
be equivalent to finding the sufficiently accurate Hermite 
quadratures. On a regular Cartesian grid, this task was 
solved systematically [20,29]. To fully recover the Galilean 
invariant Navier-Stokes momentum equation, sixth order 
accurate quadratures are required. It can be verified that 
with the D2Q17 model [17], the diagonal component of 
the sixth tensor is not isotropic, causing the Galilean in- 
variance not fully restored in the corresponding terms. 
Through exhaustive search, it can be concluded that to 
have the sixth order isotropy with velocities forming a 
regular lattice, speed-3 velocities must be included. The 
minimum velocity set found on a regular grid with sixth 
order isotropy is the 39- velocity model in three dimensions 
and the 21-velocity model in two dimensions. A different 
velocity set with sixth order isotropy has been given by 
Chikatamarla et al [26]. However, smaller velocity sets 
not entirely on a regular grid are also available [20]. 



third-order equilibrium distribution function, the diffusiv- 
ity varies with the translational velocity. Furthermore, 
this dependence of velocity also depends on the back- 
ground temperature. As shown in Fig. 2b, when fourth- 
order equilibrium distribution function is used, the ther- 
mal diffusivity does not depend on the translational ve- 
locity anymore and agrees well with the theoretical value 
for different back ground temperatures. The Galilean in- 
variance is completely recovered. However, the stability 
of the model depends one background temperature and 
translational velocity. Hence, in Fig. 2, there are no data 
for some high velocity cases. Detailed discussion of the 
instability is beyond the scope of this paper. 
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